-A223  547 


OTIC  f,le  copy 


NCEL 

Technical  Note 


_ N-1807 

February  1990 

By  L.  J.  Malvar  and  G.  E.  Warren 

Sponsored  By  Naval  Facilities 
Engineering  Command 


o 

< 


MIXED  MODE 
PROPAGATION  IN 


CRACK 

CONCRETE 


i  ! 


ff'r 

V 

ABSTRACT  Two  smeared  crack  approaches  to  fracture  of  con¬ 
crete  in  mixed  mode  are  implemented  in  two-dimensional  nonlinear 
concrete  elements:  (1)  tensile  stress  transfer  across  cracks  and 
(2)  tensile  plus  shear  stress  transfer  across  cracks.  To  corroborate 
the  analytical  model  a  notched  beam  under  mixed  mode  loading  is 
then  analyzed.  In  both  cases,  the  stiffnesses  normal  and  parallel 
the  crack  were  modified  to  insure  a  positive  definite  stiffness 
matrix.  Stresses  were  corrected  and  set  as  functions  of  the  crack 
slip  and  crack  opening.  Equilibrium  iterations  were  implemented  to 
redistribute  stress.  In  both  cases,  acceptable  agreement  was  found 
between  analytical  predictions  and  experimental  results.  The  con- 
siaeration  of  shear  stress  transfer  yielded  better  predictions,  but 
requires  consideration  of  a  non-sym metrical  stiffness  matrix,  -pi  " 
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PURPOSE 


The  Naval  Facilities  Engineering  Command  (NAVFAC)  through  the  Naval 
Civil  Engineering  Laboratory  (NCEL)  has  initiated  a  project  to  develop 
fracture  mechanics  methodology  for  design  appl  ication  of  reinforced  con¬ 
crete  elements  in  tensile  and  shear  stress  states.  In  a  preceding  study 
analytical  modeling  methodology  of  Mode  1  (opening)  was  detailed  in  two 
and  three  dimensions,  and  Mode  II  crack  propagation  (shearing)  addressed 
In  this  report  Mode  II  modeling  methodology  is  developed  and  a  benchmark 
mixed  mode  problem  is  analyzed.  This  report  supports  the  project 
"Fatigue  and  Fracture  of  Concrete"  in  the  NAVFAC  6.1  Basic  Research 
Program  YR023-03-01,  Structural  Modeling. 

The  modifications  implemented  in  the  computer  program  AD1NA  have 
been  compiled  in  the  appendixes. 

INTRODUCTION 

Although  it  is  generally  recognized  that,  crack  initiation  in  con¬ 
crete  occurs  in  Mode  I  (opening),  crack  propagation  is  more  likely  to 
take  place  in  mixed  mode,  i.e.,  involving  Mode  I  and  IT  (shearing),  or 
III  (tearing). 

Mixed  mode  crack  propagation  involves  considering  the  transfer  of 
tensile  and  shear  forces  across  cracks.  Constitutive  relations  repre¬ 
senting  the  transferred  stresses  were  evaluated  (Ref  1).  In  the  present 
report  these  constitutive  relations  are  implemented  in  a  general  purpose 
finite  element  program  developed  by  ADINA  RM)  Inc.  (Ref  2).  A  benchmark 
experiment  by  Arrea  and  Tngraffea  (Ref  3  and  4)  is  then  modeled,  with 
and  without  considering  transfer  of  shear  stresses. 


PROBLEM 

The  mixed  mode  problem  considered  is  depicted  in  Figure  1,  uni 
concrete  properties  used  are  reported  in  Table  1.  In  many  case <  the 
problem  was  approached  without  considering  shear  transfer  across  the 
crack.  Initial  attempts  at  modeling  the  shea?  transfer  using  a  constant 
shear  retention  factor  3  (typically  3  <  0.1)  yielded  results  with  almost 
no  softening  after  peak  load  (Ref  3  and  6).  Better  representations  were 
obtained  either  assuming  the  existence  of  a  Mode  II  fracture  energy 
(Ref  7),  or  using  a  predetermined  crack  path  (Ref  R). 

In  this  study  the  consideration  of  a  shear  transfer  model  is 
attempted  and  its  effects  observed. 


TENSILE  STRESS  TRANSFER 


The  transfer  of  tensile  stresses  across  a  crack  had  already  been 
implemented  with  a  smeared  crack  approach  (Ref  9)  using  the  Crack  Band 
Model  (CBM)  (see  Appendixes  A,  B,  and  C).  This  tension  softening  behav¬ 
ior  involved  a  negative  stiffness,  C  ,  for  the  cracked  element.  The  CBM 
was  implemented  assuming  zero  stiffness  (actually  a  very  small  value  was 
used  to  avoid  a  singular  stiffness  matrix)  and  then  resetting  the  stresses 
as  a  function  of  the  crack  opening.  These  stresses  are  then  redistributed 
during  equilibrium  iterations.  The  stress  transferred  versus  crack  width 
relationship  is  tabulated  in  Table  2  (Ref  1)  in  nondimens ional  form. 

The  fracture  energy,  G^,  is  related  to  the  stress  versus  displacement 
relationship  by: 

Gf  =  Wo  ft  (O0/ft  d(w/wo} 

where  w  =  crack  width  or  crack  opening 

a  =  stress  transferred  at  crack  width  w 

Wq  =  crack  width  beyond  which  no  stress  is  transferred 

ft  =  tensile  strength 

In  a  first  analysis  on  the  mixed  mode  problem  the  CBM  alone  was 
used  in  order  to  evaluate  the  importance  of  considering  shear  transfer. 

The  latest  version  of  ADINA  (Ref  2)  acknowledges  the  importance  of 
strain  softening  by  including  a  linear  stress  release  as  a  function  of 
strain  after  cracking  (Ref  10).  However,  this  stress  release  is  not 
explicitly  linked  to  fracture  energy,  and  the  authors  have  shown  that 
linearizing  the  highly  nonlinear  post  peak  stress  versus  strain  rela¬ 
tionship  negatively  affects  results  (Ref  9).  Hence,  this  feature  of  the 
latest  ADINA  version  was  not  used  in  the  present  project. 


ARC-LENGTH  PROCEDURE 

The  solution  of  the  finite  element  incremental  equations  of  motion 
was  first  attempted  using  the  spherical  arc- length  and  the  constant  incre¬ 
ment  of  external  work  procedures  described  in  Reference  11.  The  post-peak 
numerical  analysis  of  this  experiment  has  shown  to  be  highly  unstable 
(Ref  5).  The  adopted  approaches  did  not  yield  converged  equilibrium 
states  past  peak  load  and,  thus,  were  modified. 

A  type  of  indirect  displacement  control  (Ref  12)  was  then  adopted: 
in  the  arc- length  procedure  the  norm  of  displacement  (involving  all 
nodal  points)  was  replaced  by  the  distance  between  the  two  points  at  the 
edges  of  the  notch  (Appendix  C).  The  vertical  component  of  this  distance 
is  referred  to  as  CMSD  (Crack  Mouth  Sliding  Displacement).  During  the 
test,  the  CMSD  is  a  monotonically  increasing  parameter  that  stabilized 
the  algorithm.  In  the  experiment,  the  CMSD  had  been  used  as  feed-back 
control  parameter. 
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FAILURE  ENVELOPES 


The  failure  envelopes  used  in  ADINA  (Ref  13)  are  largely  based  on 
biaxial  concrete  strength  experimental  results  by  Kupfer  et  al.  (Ref  14). 
In  the  plane  stress  analytical  model,  the  crack  path  showed  sensitivity 
to  the  tensile  envelope  representation  close  to  the  tens ion/tens ion  zone 
(o^>0,  02>®»  °3=®)  (Figure  2).  The  existing  linear  envelope  in  the 
tension/compression  zone  was  then  modified  to  better  match  experimental 
results.  The  following  power  relationship  was  used: 


where  =  uniaxial  cut-off  tensile  stress 


o*  =  uniaxial  cut-off  tensile  stress  under  multiaxial  conditions 
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uniaxial  compressive  failure  stress  under  multiaxial 
conditions 

principal  stress  in  direction  i  at  time  t 


1 


if  o'c  >  8000  psi  (563  kp/cm^) 


=  1  +  0. 0002 (8000 -o'  )  if  o'  <  8000  psi 


Both  linear  and  power  envelopes  are  showjj  in  Figure  2,  together  with 
Kupfer  et  al.  results  for  Bp  =  315  kp/cm  (4450  psi).  Bp  is  the 
uniaxial  compressive  strength  of  50  by  50  by  200  mm  (2  x  2  x  7.9  in.) 
prisms.  The  current  ADINA  addresses  this  deficiency,  but  corrects  it  in 
a  different  way  (Ref  10): 


From  Figure  2  it  is  apparent  that  the  present  modification  yields  a 
better  match.  Program  modif ications  are  reported  in  Appendix  D. 

FINITE  ELEMENT  MODEL 

The  finite  element  mesh  used  is  depicted  in  Figure  1.  Loads  of 
0.13P  and  P  were  applied  at  points  A  and  B,  respectively.  In  the  com¬ 
puter  program  this  is  accomplished  using  an  automatic  step  incrementa¬ 
tion  method,  where  the  level  of  externally  applied  loads  is  adjusted 
automatically.  In  the  experiment,  a  single  total  load  of  1.13P  was 
applied  on  a  steel  beam  bearing  on  rollers  at  points  A  and  B.  The  point 
of  application  of  that  total  load  will  be  referred  as  point  C. 
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SHEAR  TRANSFER 


Cracks  in  reinforced  concrete  are  able  to  transmit  shear  forces 
across  crack  faces.  This  transfer  is  traditionally  neglected  on  the 
assumption  that  this  would  be  a  conservative  simplification.  However, 
Bazant  et  al.  showed  that  this  assumption  can  be  an  over  simplification 
(Ref  15  and  16).  Crack  dilation  occurs  with  shear  slip.  However,  crack 
dilation  is  prevented  by  forces  normal  to  the  crack  faces,  which  will 
have  to  be  compensated  by  t ensiles  forces  in  the  reinforcement  across 
the  crack. 

Shear  stresses  can  be  transferred  across  a  crack  in  three  ways: 

(1)  aggregate  interlock  as  a  result  of  the  roughness  of  the  crack  faces, 

(2)  dowel  action  or  shear  resistance  of  the  reinforcement  across  the 
crack,  and  (3)  the  axial  tensile  force  component  in  the  reinforcement 
oblique  to  the  plane  of  cracking. 

For  members  with  low  reinforcement  and  for  small  crack  widths, 
aggregate  interlock  is  the  main  mechanism  of  shear  transfer.  Tests 
carried  out  on  beams  without  web  reinforcement  showed  that  aggregate 
interlock  accounted  for  up  to  75  percenL  of  the  shear  transfer  (Ref  17). 
Hence,  most  attention  will  be  given  to  this  first  mechanism  of  transfer. 


SHEAR  TRANSFER  MODEL 


Three  accepted  empirical  models  which  represent  the.  nonlinear  rela¬ 
tionships  between  shear  stress  and  slip  are:  the  Rough  Crack  Model  (RCM) 
in  its  original  form  (Ref  11),  or  in  a  modified  form  (MRCM)  (Ref  18), 
and  the  Two-Phase  Model  (TPM)  (Ref  19  and  20).  The  constitutive  laws  of 
the.  MRCM  are  as  follows: 
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is  the  crack  stiffness  matrix. 


The  derivation  of  B  is  shown  in  Appendix  F. . 


fMPT.EMF.NTAT  I  ON  fN  FINITE  ELEMENT  PROGRAM 

Transfer  of  shear  stresses  was  implemented  by  combining  the  MRCM 
and  the  CRM.  The  incremental  flexibility  matrix  due  to  the  solid 
concrete  and  including  strain  softening  in  tension  is  given  by  (Ref  21): 
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where  p  =  Poisson's  ratio. 


(4) 


In  addition,  since  we  assume  strain  softening  in  tension  to  be  present., 
the  slope  0  of  the  strain  softening  branch  has  to  be  taken  into  account . 
The  crack  stiffness  is  then: 
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For  very  small  values  of  the  crack  opening,  C  is  large,  but 


almost  zero;  whereas,  when  the  crack 
opposite  holds. 


opening  reaches  about  0 .  1  mm , 
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The  incremental  stiffness  matrix  can  he  obtained  as  follows: 
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This  yields  an  incremental  stiffness  matrix  which  is  not  sym¬ 
metrical  and  is  not  guaranteed  to  be  definite  positive. 

Since  ADINA  considers  only  symmetric  matrices,  the  solution  was 
attempted  using  a  modified  stiffness,  and  then  correcting  the  stresses 
at  every  iteration  for  each  load  step  increment.  Tt  was  then  assumed 
that : 


V  =  0 

and  to  insure  definite  positiveness 
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where  0  is  a  small  positive  number. 


HODET.  REPRESENTATION 

To  evaluate  the  effects  of  shear  transfer,  a  100  by  100  by  100  mm 
(4x4x4  in.)  concrete  finite  element  was  first  cracked  in  tension, 
then  sheared  in  the  perpendicular  direction  (Figure  3),  in  displacement 
control.  Given  the  nodal  displacements,  strains  at  the  Gauss  points  are 
evaluated,  then  an  iterative  process  determines  crack  slip,  crack  dilata¬ 
tion,  and  concrete  deformation,  using  formulas  (1),  (21,  C,  ,  and: 

s 


The  model  behavior  is  predicted  using  all  three  formulations  (RCM, 
MRCM,  TPM).  For  each  case,  Figure  4  shows  the  shear  and  normal  loads 
transferred.  The  TPM  values  were  capped  to  the  maximum  predicted  by  the 
RCM.  It  is  observed  that  the  dilatancy  induces  vertical  compression 
(along  the  z  axis).  If  reinforcing  bars  perpendicular  to  the  cracks 
were  present,  the  dilatancy  would  increase  the  tension  in  the  bars  at 
the  crack  locations. 

From  Figure  4  it  is  apparent  that  all  three  models  yield  very  simi¬ 
lar  shear  transfer  capacity,  but  the  normal  stress  due  to  dilatation  is 
significantly  higher  for  the  RCM.  Since  more  normal  stress  experimental 
data  appears  to  back  the  MRCM  and  TPM,  the  RCM  was  discarded.  In  the 
mixed  mode  analysis,  the  more  recent  MRCM  formulation  was  chosen,  since 
it  presents  no  discontinuity  in  the  stress  gradient. 

RESULTS 

In  order  to  evaluate  the  importance  of  modeling  stress  transfer 
across  cracks,  the  analytical  model  was  first  run  with  no  transfer, 
i.e.,  assuming  total  stress  release  right  after  cracking.  Since  the 
standard  algorithms  did  not  converge,  the  indirect  displacement  method 
was  used,  with  a  very  low  fracture  energy  (0.0002  N/mm)  equivalent  of  a 
sudden  stress  release.  Results  for  this  first  run  are  shown  in  Figure  5. 

The  analysis  was  then  carried  out  considering  only  tensile  stress 
transfer  across  the  cracks  (CBM).  Finally,  the  MRCM  was  added  and  a  new 
analysis  completed  (CBM+MRCM).  Results  for  both  cases  are  shown  in  the 
form  of  load  versus  CMSD  (Figure  5),  and  load  versus  vertical  displace¬ 
ment  at  point  C  (Figure  6).  The  vertical  displacement  at  point  C  was 
derived  by  linear  interpolation  of  the  vertical  displacements  of  points 
A  and  B.  Data  points  indicating  the  reported  range  of  experimental 
results  (Ref  3)  are  shown  in  Figure  5. 

Convergence  of  the  arc- length  algorithm  was  only  obtained  for  care¬ 
fully  chosen  control  parameters.  These  parameters  control  the  size  of 
the  step  in  the  load-CMSD  space  (ALFA),  the  maximum  number  of  iterations 
allowed  for  each  time  step  (ITEMAX),  the  maximum  displacement  at  control 
point  E  (DISPP),  and  energy  convergence  criteria  (ETOL)  (Ref  11).  In 
each  case  they  were  respectively: 


Parameter 

No  Transfer 

CBM 

CMB+MRCM 

ITEMAX 

30  -4 

ETOL 

10  6 

10  6 

5.10  4 

DISPP 

-0.015 

-0.015 

-0.015 

ALFA 

0.4 

0.4 

0.5 

The  crack  pattern  for  the  last  loading  step  is  indicated  in 
Figure  7  (CBM  case).  Figure  8  shows  the  deformed  shape  obtained  for  the 
last  step  (CBM  case). 


DISCUSSION 


Figure  6  indicates  that  the  displacement  at  point  C  presents  a 
sharp  snap-back  past  peak  load.  This  explains  why  displacement  control 
at  that  point  cannot  yield  the  post  peak  response.  The  displacement  at 
both  points  A  and  B  shows  a  similar  behavior,  which  explains  why  the 
norm  of  displacement  in  the  arc-length  procedure  was  unsuccessful. 

Figure  5  shows  that  considering  tensile  stress  transfer  alone  yields 
a  conservative  behavior  prediction.  The  maximum  load  is  underestimated 
by  about  20  percent,  and  the  post  peak  load  carrying  capacity  is  lower. 
However,  the  shape  of  the  strain  softening  portion  is  similar.  A  higher 
value  of  would  yield  a  better  match  to  the  experimental  peak  load  and 
post  peak  response  (Ref  8). 

The  crack  pattern  (Figure  7)  still  differs  from  the  reported 
experimental  crack  path.  It  was,  however,  observed  that  a  small  vari¬ 
ation  in  the  mesh  size,  or  initially  larger  load  step  sizes,  would 
affect  the  path  or  result  in  bifurcation  points.  Similarly,  stiffer 
bearing  plates  would  bring  the  crack  path  closer  to  the  notch  plane. 

The  crack  path  would  easily  follow  any  of  the  different  directions 
indicated  in  Figure  9.  Th’s  would  explain  the  discrepancies  in  crack 
paths  found  by  different  authors  (Ref  22,  23,  24,  and  25)  (using  a 
similar  but  symmetrical  specimen).  For  example,  the  experimental  crack 
path  obtained  in  Reference  23  coincides  with  the  analytical  crack  pat¬ 
tern  shown  in  Figure  7. 

Should  tensile  stress  transfer  not  have  been  considered,  the. 
maximum  load  carrying  capacity  of  the  analytical  model  would  have  been 
reached  as  soon  as  the  first  tensile  cracks  formed  (around  50  kips) 
(Figure  5).  This  is  obviously  an  inadequate  representation  of  t^e 
experimental  behavior. 

Transfer  of  both  tensile  and  shear  stress  is  considered  best  in 
matching  experimental  behavior.  The  peak  load  is  higher  and  the  post 
peak  behavior  is  closer  to  experimental  results.  However,  in  order  to 
obtain  the  complete  post  peak  behavior,  a  nonsymmetr i ca 1  stiffness 
matrix  would  have  to  be  considered.  This  would  present  additional 
difficulties,  such  as  (1)  implementation  in  a  new  program  with  a 
nonsymmotrical  solver,  and  (2)  increase  in  computation  time.  The 
increased  accuracy  has  to  be  weighed  against  the  increased  cost  in 
Implementing  shear  transfer.  In  this  case,  the  crack  pattern  remained 
similar  to  the  previous  one. 


CONCLUSIONS 

The  consideration  of  shear  stress  transfer  across  the  propagating 
cracks  yielded  a  better  prediction  of  the  experimental  results.  However, 
the  resultant  stiffness  matrix  is  nonsymmetr ical  and  would  require  imple¬ 
mentation  in  a  program  with  a  nonsymmetric  solver.  This  would  enhance 
the  convergence  of  the  indirect  displacement  control  '■lgorithm. 

The  exclusive  consideration  of  tensile  stress  transfer  yielded  good 
results  up  to  peak  load.  Beyond  tin’s  point,  the  loads  are  underestimated, 
although  the  shape  of  the  unloading  branch  matches  the  experimental  trend. 
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Tli  i s  could  he  an  acceptable  represent at  ion  of  mixed  mode  behavior  as 
long  as  it  is  kept  in  mind  that  a  conservative  post  peak  behavior  will 
be  obtained.  Finally,  it  was  shown  that  inadmissible  results  are  obtained 
if  both  tensile  anil  shear  stresses  are  assumed  to  completely  vanish  upon 
cracking. 
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1  1 


Tabic  1.  Concrete  Properties 


Fracture  Energy 
Compressive  Strength 
Tensile  Strength 
Modulus  of  elasticity 


=  0.055  N/mm 

f'c  =  45.5  N/mm2 
f  =  2.80  N/mm 

F.  =24.8  GPa 


Table  2.  Stress  -  Crack  Width  Relationship 


w/w 

o 

a/^ 

0.00 

1.0000 

0.05 

0.7082 

0.  10 

0.5108 

0.  15 

0.3817 

0.20 

0.2986 

0.25 

0.2446 

0.  30 

0.2080 

0.40 

0.  1596 

0.60 

0.0904 

0.80 

0.0361 

1 . 00 

0.0000 

Tension-compression  failure  envelope 
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ROUGH  CRACK  MODEL 


Y  REACTION  (kN) 


~ Two-Phase  Model 
~+~  Rough  Crack  Model 
Modified  Rough  CM 


u  u.  u6 


Z  DISPLACEMENT 
Y  DlSPL.  (mm) 


(a)  Y  Reaction  (kN) 


Z  REACTION  (kN) 


(b)  Z  Reaction  (kN) 


Figure  4.  Shear  transfer  models. 


Z  DISPLACEMENT 
Y  DlSPL.  (mm) 


Figure  5.  Load  versus  CMSD  plots. 


Figure  7.  Crack  pattern. 


Figure  9.  Alternate  crack  paths 


Appendix  A 


CBM,  2-i) 


The  following  are  changes  implemenLed  in  the  two-dimensional  element 
formulation;  namely,  in  the  subprograms  TODMFF..F77  and  ELT2D4.F77  for 
tensile  stress  transfer.  Although  not  mentioned  in  this  report,  formula¬ 
tion  for  the  Rotating  Crack  Model  is  also  included  (see.  Ref  1). 

Version  84.NL3  of  the  computer  program  ADTNA  was  used.  Text  record 
identifiers  were  inserted  in  Columns  72  to  80,  using  HISTORIAN  PLUS 
Software  Control  System. 


CHANGES  IN  TODMFE . F77 

Changre  at 
or  after i 


1  IDWAS/  0,  0,  0,18,18,  0,10,15,15,5?, 55,  0,  0,26,6*0/,  T00MFE9? 

COMMON  /SOFT/  ISCODE ,WWCC ,ELWW ,GGFF ,DDAA , IRCM  TDFE  4Z 

IF  ( MODEL . EQ . 5  1  READ! I IN , 1005  I  ISCOOE ,WWCC ,ELWW ,GGFF ,DDAA , IRCM  TDFE  101 

1005  FORMAT  (I5.4F10.0I  TDFE 121 9 

COMMON  /SOFT/  ISCODE ,WWCC ,ELWW ,GGEF ,DDAA , IRCM  MATRT214 

WRITE  (6,22561  ISCODE ,WWCC ,ELWW ,GGFF ,DDAA  MATRT244 

IF  (IRCM. EQ. 01  WRITE  (6,22571 
IF  (IRCM.GE.il  WRITE  (6,22581 

2256  FORMAT* /58H  ( BB  I  CODE  FOR  TENSILE  STRESS  TRANSFER  ,T5.  MATRT726 


1  /58H  1 -LINEAR  SOFTENING 

2  /58H  2 “CORNEL ISSEN ' S  SOFTENING 

5  /58H  SOFT  BAND  WIDTH  (WWCCI  ,F10.5. 

4  /58H  SOFT  ELEMENT  WIDTH  (ELWWI  ,FI0.5, 

5  /58H  FRACTURE  ENERGY  ( GGFF  1  .riO.O, 

6  /58H  MAXIMUM  AGGREGATE  SIZE  (DDAAI  ,ri0.5l 

2257  FORMAT! /41H  ONLY  PERPENDICULAR  CRACKS  ALLOWED  I 

2258  FORMAT! /41H  ROTATING  CRACK  MODEL  IS  USFD  I 


A-l 


CHANGES  IN  ELT2D4 .  F  77 

IDW=18*ITW0 

DIMENSION  PROP!  1  I.WA!  18,1  I  ,YZ<  1  l,NOD5(  1  l.NODS!  1  1.TEMPV1I  1  t 

DO  10  1=1,18 

COMMON  /SOFT/  ISCODE ,WWCC ,ELWW ,GGFF ,DDAA , IRCM 

1  CRKSTR!  6  I, STRESS!  4  I, STRAIN!  4  I, C(  4, 4  1  ,NOOS<  1  I  ,TEMPV1(  1  I 

2  TEMPV2!  1  t  ,YZ(  1  I  ,N005<  1  I  ,WA(  1  >  .DUMWAI  18  » 

DO  1  1=1,18 

IF  (IRCM.GE.l  .AND.  ANGLE . LT . 5. 61D2 t  GOTO  II 

GO  TO  14 
IS  CONTINUE 

CALL  CRAKID  ( STRESS .STRAIN ,PGRAV .CRKSTR ,RKLD ,RKUN ,GU) , SPSS , 

1  ANGLE, EP, NUMCRK, MODEL, 1  I 

14  CONTINUE 

47  CALL  DCRACK  I C ,SIG .ANGLE .MODEL , ITYP2D .NUMCRK , 1 , 1 ,CRKSTR ! 

CALL  DCRACK  ( C .STRESS ,ANG ,MOOEL , ITYP2D .NUMCRK ,1,2 .CRKSTR  I 

CALL  DCRACK  ( C .STRESS .ANGLE .MODEL , ITYP2D .NUMCRK , 2 , 2 .CRKSTR  I 

CALL  DCRACK  ( C .STRESS .ANGLE .MODEL , ITYP2D .NUMCRK , 1  , 2 .CRKSTR * 

CRKSTR!  4  l=EP<  1  I 
CRKSTR!  5  l=EPI  2  > 

CRKSTR!  6  l=EPt  5  ! 

CALI.  DCRACK  I  C  .STRESS  .ANGLE  .MODEL ,  ITYP2D  .NUMCRK  ,  1  , 2  .CRKSTR  I 
CALL  DCRACK  ! C .STRESS .ANGPRI .MODEL , ITYP2D .NUMCRK , 1 ,2 .CRKSTR  I 
CALL  DCRACK  ( C .STRESS ,ANG .MODEL , ITYP2D , NUMCRK .2 . 1 .CRKSTR  I 
DO  210  1  =  1 ,18 

COMMON  /SOFT/  ISCODE  ,WWCC  .EI.WW  ,GGFF  ,DDAA  ,  I RCM 

DIMENSION  STR!  4  I  .EPS!  4  I  .CRKSTR!  6  I  ,SF1  (  1  )  ,SPS1  (  1  •  .SPS7I  1  I  .SPSS!  1 
IF  ! IRCM.GE.l  I  GO  TO  11 

GO  TO  107 


A -2 


ELT2D458 

ICDM0D16 

ICDM0D26 

CDMOD  50 

,  CDMOD  51 

CDMOD  54 

CDMOD  66 

CDMOD 1 35 

CDMOD 150 

CDMOD 2 70 

CDMOD 302 

CDMOD 350 

CDMOD 3 74 

CDM0D415 

CDM00422 

CDMOD427 

CDM00590 

CDM0D596 
CRAKID1 3 

,  CRAKID15 

CRAKID16 

CRAKID45 


c 

1 1  IF  ( KKK . GE . 2  1  GO  TO  12 
C 

C  FIND  DIRECTION  OF  PRINCIPAL  STRAINS 
C 

AA=<  EPS<  1  I  4  EPS!  2  I  1*0.5 

BB=(  EPS(  1  I  -  EPS(  2  11*0.5 

CC=SQRTI  BB*BB  4  EPSI  3  >*EPS(  3  1  1 

EPSLI  1  1=AA  4  CC 

EPSL<  2  l=AA  -  CC 

EPSLI  3  1=0. DO 

EPSLI  4  l=EPS(  4  ) 

ANGLE=4.5D1 

IF  (EPSI  3  I.EQ.O.DO  I  ANGLE=0.1D-3 
IF  (ABSIBB1.LT.0.1D-6I  GO  TO  12 
DUM-ABSI  EPSI  3  l/BB  1 
ANGLE=57. 296*ATANI DUM  1 
C 

IF  (BB.LT.O.DO  .AND.  EPSI 3  I . GT . 0 . DO  1  ANGLE  =  180 .  -  ANGLE 

IF  (BB.LT.O.DO  . AND .  EPSI 3  I . LE . 0 .DO  I  ANGLE  =  180 .  4  ANGLE 

IF  (BB.GT.O.DO  .AND.  EPSI 3  I . LE . 0 .DO  1  ANGLE=360 .  -  ANGLE 

ANGLE=ANGLE/2. 

C 

C  FIND  STRESSES  PERPENDICULAR  AND  PARALLF.I.  TO  CRACK 
C 

12  CONTINUE 
PI=4.D0*ATAN( 1 .DO  1 
TANG=ANGLE 

IF  (TANG.LT.-5.41D2  I  TANG=TANG  4  722. 

IT  (TANG  .LT.  (-1.8D2I1  TANG=TANG  4  361. 

IF  I  TANG . GT . 1 .802  I  TANG=TANG  -  180. 

GAM  =  2  .  *ABS(  TANG  l«PI/180. 

SG=SINI GAM  1  CG=COS( GAM  I 
IF  I KKK . EQ . 3  1  GO  TO  107 
C 

R1 1  =  1  STRI  1  I  4  STRI  2  11*0.5 
R1 2  =  1  STRI  1  1  -  STRI  2  1  1*0 . 5 
SIGPI  1  I=R1 1  4  R12«CG  4  STRI  3  1*SG 
STGPI  2  I=R1 1  -  R1 2*CG  -  STRI  3  1*SG 
SIGPI 3  1=0. DO 
SIGPI  4  l=5TRI  4  1 
C 

IF  (KKK.EQ.2I  RETURN 

COMMON  /SOFT/  ISCODE  ,WWCC  ,FLWW  ,GGFF  >DDAA ,  IRCT1 
DIMENSION  Cl  4  ,4  I  ,SIG(  4  I  ,DI  4  ,4  I  ,T(  4  ,4  1  ,DSTGI  4  1  .CRKSTRI  6  I 


DCRACK  8 

DCRACK  9 


A-3 


IF  ( IRCM.EQ.O  I  GO  TO  12 
IF  (  EPI  1  I.NE.EPI2I  I 


DCRACK59 


1  C(  3  >  5  I  =  (SIGPI  1  l-SIGPI  2  I  I/I  2*(EPI  1  I  EPI  2  I  I  I 
IF  (  EP(  I  I .  EQ  .EP(  2  I  I  Cl  3 , 3  I  =  1.0-8 
12  CONTINUE 

C  RELEASE  APPROPRIATE  STRESSES 

C 

98  NF=NUMCRK  ♦  1 

GO  TO  (140,120,110,155,100,100,1001,  NF 
100  CALL  DSOF  ( 4 ,SIGP ,FALSTR ,EP ,CRKSTR  ,E ,VNU ,SIGMAT ,SIGMAC  I 
IF  (NUMCRK  -  51  140,120,110 

110  CALL  DSOF  ( 2  ,SIGP  ,FALSTR ,EP ,CRKSTR ,E , VNU  >SIGMAT , SIGH AC  I 
120  SIGP(  3  )=SISP(  3  I 

CALL  DSOF  ( 1 ,SIGP ,FALSTR ,EP .CRKSTR ,E ,VNU ,SIGMAT ,SIGHAC  I 
C 

C  ROTATE  STRESSES  TO  GLOBAL  AXES 

SUBROUTINE  DSOF  I IJ ,SIGP .FALSTR ,EP ,CRK5TR ,E , VNU ,SIGMAT ,SIGMAC  I 
IMPLICIT  DOUBLE  PRECISION  (  A-H.O-Z  I 
COMMON  /SOFT/  ISCODE ,WWCC ,ELWW ,GGFF ,DDAA , IRCM 
DIMENSION  SIGPI4  I  ,EP(  4  I  .CRKSTRI  6  I  ,CORN(  11,31 
IF  (CRKSTRI  IJ  I.GT. 0. DO  I  GOTO  5 
SIGPI  IJ  l=FALSTR 
RETURN 
5  CONTINUE 
C 

DATA  ( CORN! 1,1  1 , 1  =  1 , 1 1  I/O . , . 05 , . 1  , . 15 , . 2 , . 25 , . 3 , .4 , . 6 , . 8 , 1 . 0/ 
DATA  (CORN! 1,2  1,1  =  1 ,11  1/1 . 7082 , . 51 08 , . 381 7, . 2986 , . 2446 , 

1  .2080, .1596, .0904, .0361, 0.0/ 

JJ=IJ 

IF  (JJ.EQ.4I  JJ=3 
KK=JJ*3 
EEPP=EP( IJ  I 

IF  (EPI  IJ  I.GT.  CRKSTRI  KK  I  I  CRKSTRI  KK  1=EP<  IJ  I 
IF  (EPI  IJ  I.LT.  CRKSTRI  KK  I  I  EEPP=CRK5TR(  KK  I 
ISS=ISCODE-Z 
IF  I ISS  I  10,20,30 
C 

10  CONTINUE 

EETT=l/( l/E-l 2*GGFF  I/I SIGMAT«*2*WWCC I  I 

SIGPI  IJ  1=FALSTR»EETT*(  EEPP-CRKSTRI  JJ  I  I 

IF  (EPI  I J  I.LT.  CRKSTRI  KK  I  I  SIGPI  IJ  l=EPI  IJ  l/FEPP»SIGPI  TJ  I 

IF  (SIGPI  IJ  I.GT. FALSTR  I  SIGPI  IJ  l=FALSTR 

IF  (SIGPI  IJ  I.LT. 0. DO  I  SIGPI  IJ  1=0. DO 

SIGPI  3  1=0. DO 

RETURN 


A-4 


DCRAC204 
DCRAC205 
DCRAC206 
DCRAC207 
DCRAC208 
DCRAC209 
DCRAC210 
DCRAC21 1 
DCRAC212 
DCRAC21 3 
DCRAC214 

CDMOD620 


20  CONTINUE 

EO=GGFF/IWWCC»0. 19704*SIGMAT I 
DO  21  1=1,11 

CORN!  1,3  !=CORN<  1,1  HCORNI  1,2  l»CRKSTR<  JJ  l/EO 
IF  ( EEPP/EO.LT. CORN! 1,5  It  GO  TO  22 

21  CONTINUE 

22  AA=<  CORN!  1-1,2  )-CORN(  1,2  I  I/I  CORN!  I  1  ,3  I -CORN!  1,3  11 
BB=CORN(  1-1,2  l-AA»CORNI  1-1,3) 

SIGPI IJ  l=FALSTR«< AA*EEPP/EO+BB 1 

IF  ( EP( IJ  ) . LT.CRKSTRI KK  I  )  SIGPI IJ  l=EPI IJ  l/EEPP»SIGP< IJ  I 
IF  (SIGPI  IJ  I.GT.FALSTR  I  SIGPI  IJ  )=FALSTR 
IF  I  SIGPI  IJ  I . LT . 0  .DO  I  SIGPI  IJ  1=0. DO 
SIGPI  3  >=0  .DO 
RETURN 
C 

30  CONTINUE 

RETURN 
C 

END 
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Append i x  B 
CBM,  3-D 


Three-dimensional  element  formulation. 


CHANGES  IN  THREDM.F77 


1  IDWAS  /  0,  0,  0,  25,25,  0 , 14 , 21  , 21  ,47 ,47 , 38 ,B«0/ , 

COMMON  /SOFT/  ISCODE ,WWCC,ELWW ,GGFF ,0DAA , 1 RCM 
IF  (MODEL. EQ. 51  READ!  UN ,  1009  )  ISCOOE  ,WMCC  ,ELWW  ,GGFF  ,DDAA ,  1RCM 
1009  FORMAT  (I5,4F10.0I 


COMMON  /SOFT/  ISCODE ,WWCC ,ELWW ,GGFF  ,DDAA , 1 RCM 


WRITE  (6,22391 


2239  F0RMAT(/38H  (BBI 

1  /38H 

2  /38H 

3  /38H 

4  /38H 

5  /38H 

6  /38H 


CODE  FOR  TENSILE  STRESS  TRANSFER, 15, 

1 =LINEAR  SOFTENING  , 

2=C0RNELISSEN‘S  SOFTENING  , 

SOFT  BAND  WIDTH  ( WWCC  I  ,F10.5, 

SOFT  ELEMENT  WIDTH  (ELWW1  ,F10.5, 

FRACTURE  ENERGY  I GGFF  1  ,F10.8, 

MAXIMUM  AGGREGATE  SIZE  I DD  A  A  I  .F10.51 


CHANGES  IN  ELT3D4.F77 

IDW=25*ITWO 


DIMENSION  PROP( 1  I  ,WA( 25,11  ,XYZ( 1  I , NOD 9 ( 1  I  ,NOOSI 1  1  ,TEMPV1 1  1  1 


DO  10  1=1,25 


1  CRKSTRI  6  I  .STRESSI  6  1  .STRAIN!  6  1  ,C(  6,6  1  ,RLMN(  3,3  I, NODS!  1  1 

1  TEMPVK  1  I  ,TEMPV2(  1  l,XYZ(  1  1  .N009I  1  1  ,WA(  1  I.DUMWAI  25  I 


Change  at 
or  after, 

THRED100 

THDFE  46 

THDFE102 

THDFU90 

MATWRT14 

MATWR243 

MATWR596 


ELT3D444 

1CM0031 6 

10100326 

,  CM 00  3D  54 

CM003055 


B-l 


DO  1  1  =  1  ,25 


CM0D3D67 


47  CALL  DCRAK 3  ( C ,SIG .RLMN .MODEL .NUMCRK , 1 , 1 .CRKSTR  I 

CALL  DCRAK 3  ( C .STRESS ,RLMN .MODEL .NUMCRK , 1 ,2 .CRKSTR  1 

CALL  DCRAK3  ( C .STRESS .RLMN .MODEL .NUMCRK , 2 ,2 .CRKSTR > 

CRKSTR!  4  l=EP(  1  ) 

CRKSTR! 5  t  =EP < 2  I 
CRKSTR(6>=EP(3» 

CALL  DCRAK3  ( C .STRESS ,RLMN .MODEL .NUMCRK , 1 ,2 .CRKSTR  I 
159  CALL  DCRAK 3  ( C .STRESS ,RLMN .MOOEL .NUMCRK , 1  , 2  .CRKSTR  I 
CALL  DCRAK3  ( C .STRESS , RLMN .MODEL .NUMCRK , 1 , 2 .CRKSTR  I 
130  CALL  DCRAK3  < C ,SIG , RLMN .MODEL .NUMCRK , 2 , 1 .CRKSTR  I 
DO  210  1=1 ,25 

DIMENSION  STR!  4  I  .EPS!  4  I  .CRKSTR!  6  I  ,SP1  (  1  I.SP5K  1  >,SP32I  1  » ,SP33l  1  I 

DIMENSION  Cl  4 .4  I ,SIG(  4  1 ,0<  4 ,4  )  ,T( 4 ,4  I  .DSIGl 4  > .CRKSTRI A  I 

C  RELEASE  APPROPRIATE  STRESSES 
C 

NF  =  IK  ♦  1 

GO  TO  (140,120,110,100,1551,  NF 
100  CALL  DSOF  3  (  5  ,SIGP  .FALSTR  ,EP  .CRKSTR  ,E  ,  VNU  .SIGMAT  .SIGT1AC  I 
110  SIGP!  6  1=SIGP(  6  I 

CALL  DSOF 3  I  2 ,SIGP .FALSTR ,EP .CRKSTR ,E , VNU .SIGMAT  .SIGMAC  I 
120  SIGP!  5  l=SIGP(  5  I 
SIGP!  4  )=SIGP(  4  I 

CALL  DSOF?  ( 1 .SIGP, FALSTR, EP, CRKSTR, E, VNU, SIGMAT, SIGMAC  ! 

C 

C  ROTATE  STRESSES  TO  GLOBAL  AXES 

SUBROUTINE  DSOF3  I  I J .SIGP .FALSTR ,EP .CRKSTR ,E .VNU .SIGMAT .SIGMAC  I 
IMPLICIT  DOUBLE  PRECISION  I  A-H.O  Z  I 
COMMON  /SOFT/  ISCODE ,WWCC ,ELWW ,GGFF .DDAA.1RCM 
DIMENSION  SIGP!  4  !  ,EPI  4  I  .CRKSTR!  6  I  .CORN!  11,31 
C 

IF  (CRKSTR! IJ  I.GT.O.DO  I  GOTO  5 

SIGP!  IJ  >=FALSTR 

RETURN 
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CMOD3261 

01003206 

CMOD3340 

01003362 


01003363 

01003414 

CMOD3420 

CMOD3561 

O10D3567 

,  CRAKID1 5 

DCRACK  9 

DCRAK 165 
DCRAK 166 
DCRAK167 
DCRAK 168 
DCRAK 169 
DCRAK 170 
DCRAK 171 
DCRAK 172 
DCRAK 173 
DCRAK  174 
DCRAK 175 
DCRAK 176 

010D3590 


5  CONTINUE 

DATA  ( CORN! 1 , 1 I , 1  =  1 , 1 1  I/O . , . 05 , . 1 , . 15, . 2 , . 25 , . 5 , . 4 , . 6 , . 8 , 1 . 0/ 
DATA  (CORN) 1 ,2  1 ,1=1 ,1 1  1/1 . , . 7082, . 51 08 , . 581 7 , . 2986, . 2496, . 2080, 
1  . 1596, .0904, .0561 ,0.0/ 

JJ=IJ 
KK=JJ»5 
EEPP=EP(  IJ  1 

IF  IEPI  IJ  I.GT.CRKSTRIKK  I  I  CRKSTRI  KK  l=EPI  IJ  I 
IF  IEPI  I J  I.  LT.  CRKSTRI  KK  I  I  EEPP=CRKSTR(  KK  I 

ISS=ISCOOE-2 
IF  ( ISS  I  10,20,50 

10  CONTINUE 

EE7T=l/< l/E-l 2«GGFF  l/< SIGMAT*»2*WWCC  1 1 

SIGPI  IJ  1=FALS>TR+EETT»(  EEPP -CRKSTRI  JJ  I  I 

IF  <EP<  IJ  I  .LT. CRKSTRI  KK  I  I  SIGPI  IJ  l=EPi  I J  1/EEPP*SIGPI  I J  I 

IF  ( SIGPI IJ  I . GT. FAL5TR  I  SIGPI IJ >=FALSTR 

IF  (SIGPI  IJ  I.LT.O.DO  I  SIGPI  IJ  1=0. DO 

IF  I IJ-2  I  12,11,11 

11  SIGPI  6  1=0. DO 

12  SIGPI  5  1=0  .DO 
SIGPI  4  1=0. DO 
RETURN 

20  CONTINUE 

EO=GGFF/( WWCC*0 . 1 9704*SIGMAT  I 
DO  23  1=1,11 

CORN!  1 , 5  l=CORH(  1 , 1  l*CORNI  1 , 2  l«CRKSTRI  JJ  l/EO 
IF  I EEPP/EO . LT . CORN! 1,511  GO  TO  24 
25  CONTINUE 

24  AA=I CORN!  1-1,2  l-CORNI  1,2  I  I/I  CORNI  1-1  ,5  I  -CORNI 1 ,5  I  I 
BB=CORN< 1-1 ,2  l-AA*CORN( 1-1,5  1 
SIGPI IJ  l=FALSTR»l AA*EEPP/EO*B8  I 
IF  I  SIGPI IJ  I . GT. FAL5TR  1  SIGPI  IJ  l=FALSTR 
IF  ISIGPI  IJ  I.LT.O.DO  I  SIGPI  IJ  1=0  DO 
IF  I IJ-2  I  22,21,21 

21  SIGPI  6  1=0. DO 

22  SIGPI 5  1=0. DO 
SIGPI  4  1=0. DO 
RETURN 

50  CONTINUE 
RETURN 

END 
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Appendix  C 


GENERAL  MODIFICATIONS 


The  following  are  changes  implemented  in  the  rest  of  the  program, 
namely  in  the  subprograms  ADINA.F77,  ADINI.F77  and  ADINA2.F77.  Only  the 
modified  spherical  constant  arc- length  scheme  is  allowed  and  only  Full 
Newton  iterations  without  line  search  are  carried  out.  Tf  NODQL  is 
chosen  between  3  and  100,  a  subset  of  NODQL  nodes  is  used  in  the  norm  of 
displacement.  If  NODQL  is  2,  the  distance  between  two  points  is  used 
instead  of  the  norm  of  displacement. 


CHANGES  IN  ADINA . F77 

COMMON  /DICS/  DISPM4 ,ADNOM ,ADMAX, ADCOM ,N0DQ,NDID .NEDPM4 ,N1CRL0,IARADINA189 
I  , NODQL ,NEDPML I  100,71 

IF  (IRSM4.EQ.2l  KSTOP=l  ADINA994 

CHANGES  IN  ADINI.F77 


1  , NODQL, NEDPML)  100,7 1  ADINI  33 

IRE AD  ( I IN , 1 004  I  N0DQ ,NDID .DISPM4 , ADN0M , ADMAX , ICOMA , I AR .NODQL  ADINI166 

IF  I NODQL . EQ . 0  .OR.  METH00.NE.4I  GO  TO  70 
NN0DQL=INT( FLOAT! NODQL-1  1/10.  1*1 
DO  69  1=1 ,NN0DQL 

69  READ  (IIN,  10071  (  NEDPMLI  10*<  I  - 1  UK  ,  1  t  ,K  = 1  , 1 0  I 

70  CONTINUE 

IF  (METHOD. EQ. 4 (  NEWREF=1  ADINI718 

IF  ( NODQL. EQ.O  I  GO  TO  46  ADINI801 

WRITE  (6,20671  NODQL 
DO  446  1=1 ,NNOOQL 

446  WRITE  ( 6 , 2068  •  INEDPMLf  10»l  1-1  UK,1  1.K-]  ,10  • 


1004  FORMAT  ( 215 ,3F10 . 0 ,315  I  ADIN1218 

1007  FORMAT  ( 10X5  I  ADIN1220 

2067  FORMAT  I/5X,  ADIN1577 

1 55HN0DE  SUBSET  FOR  DISPLACEMENT  NORM,  TOTAL  NODES  (NODQL >=,15,/ I 

2068  FORMAT  ( 15X.10I4  I 
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CHANGES  IN  ADINA2.F77 

1  , NOOQL, NEOPML! 100,7  I 

IF  ( NOOQL . EQ . 0  I  GO  TO  13 
DO  13  11=1, NOOQL 

NIDL=N5  -  1  ♦  I  I  NEDPML!  11,1  1-1  l*NDOF 
DO  22  IN=1 ,6 

IF  I  IDOF!  INI.EQ.O  I  GO  TO  22 
NIDL=NIDL+1 

NEDPMLI  II, IN* I  )  =  IAINIDLI 
22  CONTINUE 

13  CONTINUE 

IF( NOOQL. EQ. 01  GO  TO  14 
WRITE  ( 6 , 5000  ) 

DO  15  1=1 .NOOQL 

15  WRITE  16,5001)  NEDPML!  1,1  ) ,  (  NEDPMLI I  >K  I  ,K=2 , 7  I 

14  CONTINUE 

5000  FORMAT! /34H  NODE  SUBSET  FOR  DISPLACEMENT  NORM, 

1  /34H  NODE  EQUATION  NUMBERS  ) 

5001  FORMAT! 6X, 14, 4X, 614  ) 

C  IF  <PEOLD.GT.BIG*PEINITI  GO  TO  210 

GO  TO  230 


1  , NOOQL, NEDPML!  100,7  I 


1 


•  NOOQL, NEDPML!  100,7  1 


DUALL= 3 . *DUALL 

IF  ( NOOQL . NE . 0  I  GO  TO  500 

COMMON  /DICS/  DISPM4 ,  ADNOM  ,ADMAX  ,  ADCOM  ,NODQ,NDID  .NEDPM4  .N1CR1.0 
1  , NOOQL, NEDPML! 100,71 

IF  ( NOOQL . EQ . 2  )  GO  TO  160 
IF  I NOOQL . NE . 0  )  GO  TO  150 

150  PD=0 .DO 

DO  151  1=1, NOOQL 
DO  151  J=2 ,7 

IF  (NEDPML!  I,J  I.EQ.O  I  GO  TO  151 
PD=PD+AA(  NEDPMLI  I ,  J  I  l*BB(  NEDPMLI  I ,  J  II 

151  CONTINUE 

PD=PO»(  NALLEQ/I  NOOQL*  1  I  I 

PR=0 .DO 

RETURN 


LOADMS55 

LO ADM 103 


LOADM109 

LO ADM 268 

EQUIT254 

EQUIT259 

ASTIM423 

ASTCHE71 

ASTCH151 

ASTCH186 

IARD0PRFM14 

D0PRFM21 

DOPRFM46 
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160  PD=0 . DO 


DO  161  J=2,7 

IF  (NEDPMLI  1  ,J  I.EQ.O  I  GO  TO  161 
PD=PD  +  (AA(  NEDPMLI  1  ,J  I  l-AAl  NEDPMLI  2 ,  J  I  I  I 
1  *<  BB(  NEDPMLI  1 ,  J  I  l-BB(  NEDPHLI  2 ,  J  II  I 

161  CONTINUE 

PD=PD»<  NALLEQ/3  ( 

PR=0 .DO 
RETURN 


1 

.NOOQL 

.NEDPMLI  100,7  1 

ALSTEP  9 

1 

,  NOOQL 

.NEDPMLI  100,7  1 

ALSET  10 

NOOQL =0 

ALSET  37 

DO  2  1=1,100 

DO  2  J=1 ,7 

ALSET  39 

2  NEDPMLI  I, J  1=0 


1 


.NOOQL,  NEDPMLI  100,7  1 


NEWDAV36 
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Appendix  D 


FAILURE  ENVELOPE 

Concrete  Strength 


8000 

psi 

(563 

kp/cm2 ) 

(55.2 

MPa) 

6000 

ps  i 

(422 

kp/cm2 ) 

(41.4 

MPa) 

4000 

psi 

(281 

kp/cm2 ) 

(27.6 

MPa) 

2000 

psi 

(141 

kp/cm2 ) 

(27.6 

MPa) 

CHANGES  IN  ELT204.F77 

FALSTR=SIGMAT*( 1 .  -  ( P3/SIGMAC  l**POWR  I 

IF  (P1.GE.0.D0I  FALSTR=SIGMAT*( 1 .  -  ( PZ/SIGCP t»»POWR  1 
1  *1  1 .  -  (  P3/5IGCP  l**POWR  I 

IF  C  MODEL . EQ . 5  .AND.  SIG( Z  1 . LT . 0 .DO  I 
1  FALSTR=SIGMAT*( 1 .  -  ( SIG( Z  l/SIGMAC  l**POWR  ) 

CHANGES  IN  ELT3D4.F77 

1Z0  FALSTR=SIGMAT*( 1 .  -  < P3/SIGMAC >**POWR t 

IF  (P1.GE.0.D0)  FALSTR=SIGMAT*I  1 .  -  ( P2/SIGCP  l»*POWR  1 
1  *11.  -  < P3/SIGCP  l«»POWR  1 

IF  (MODEL. ED. 5  .AND.  SIG( 3  I . LT. 0 .DO  1 
1  FALSTR=SIGMAT*( 1 .  -  ( SIG<  3  l/SIGMAC  l*»POWR  1 


D-l 


Power 

1.00 

1.40  ' 

; 

1.80 

i 

2.20  j 

] 

PRNCPL86 

PRNCP1Z1 

CRAKI183 

CRAKI184 

PRNCPZ3Z 

PRNCP266 

CRAK31Z0 

CRAK31Z1 


Append  i-X  E 

DERIVATION  OF  CRACK  STIFFNESS  MATRIX  B 


The  MRCM  formulation  can  be  rewritten  as 


°nt  =  V'-VW'"/*1 


where: 


f  =  an  +  lr  I 


g  =  1  +  a4  r 


h  =  (Hr2)0'25 


and  by  derivation: 


f  = 


f,  = 


= 


8.  = 


h  = 


h.  = 


-3«4|6t3/«n4| 


3  f  /  3n 

3f/3t  =  3»4«t|6t/6nJ| 
3g/3n 
3g/3t 
3h/3n 
3h/3t 


-4a. ( 6  4/6  5) 
4  t  n 


4a  (6^/6  4) 
4  t  n 


2  -0  7S  n  3 

(1  +  r  )  ^(-26^  /  6  1/4 

t  n 


2  -0  7S  2 

(1+r  )  u‘ 7 ’(26/6 
t  n 


The  crack  stiffness  terms  are  then: 

Bnn  =  -ai2((-h6t/6n2-hnr)V^°nt/h2  +  r6n‘° '  5°nt/2h  rV^Btn/h) 
Bnt  =  -al2((‘htr+h/6n)V^0nt/h2  +  r'/*n\t/h) 

Btt  =  f/6n8  +  r^ft8-f8t)/«2 

Btn  =  To(-fr/(gV^)H1-\^X)[(fn8-f8n)r/g2-fV(86n2)]) 
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